1. Forecasting Lake HURON with ARMA (Direct fit)

• Recall that Lake Huron data was fit three different way (c.f. Lec 14), each yielding adequate model.

If we treat it as “stationary” series, we can try to fit ARMA(p,q) series.

  library(forecast)

## Warning in as.POSIXlt.POSIXct(Sys.time()): unknown timezone 'zone/tz/2018c.
## 1.0/zoneinfo/America/New_York'

  source("http://gozips.uakron.edu/~nmimoto/477/TS_R-90.txt")
  
  D <- read.csv("http://gozips.uakron.edu/~nmimoto/pages/datasets/lake.csv")
  Lake  <- ts(D, start=1875, freq=1)
  length(Lake)

## [1] 97

  plot(Lake, type="o")

  layout(matrix(c(1,2), 1, 2))
  acf(Lake)
  pacf(Lake)

  layout(1)

  Stationarity.tests(Lake)

## 
##     'tseries' version: 0.10-42
## 
##     'tseries' is a package for time series analysis and
##     computational finance.
## 
##     See 'library(help="tseries")' for details.

## Warning in kpss.test(A): p-value smaller than printed p-value

##        KPSS  ADF    PP
## p-val: 0.01 0.24 0.025

  Fit1 <- auto.arima(Lake, d=0, stepwise=FALSE)      # If you force ARMA(p,q) fit
  Fit1    

## Series: Lake 
## ARIMA(1,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1     ma1    mean
##       0.7665  0.3393  9.1290
## s.e.  0.0773  0.1123  0.3861
## 
## sigma^2 estimated as 0.4784:  log likelihood=-101.09
## AIC=210.18   AICc=210.62   BIC=220.48

  Randomness.tests(Fit1$residuals)

##   B-L test H0: the sereis is uncorrelated
##   M-L test H0: the square of the sereis is uncorrelated
##   J-B test H0: the sereis came from Normal distribution
##   SD         : Standard Deviation of the series

##       BL15  BL20  BL25  ML15  ML20   JB    SD
## [1,] 0.963 0.952 0.934 0.567 0.641 0.89 0.684

Analysis 1 (direct fit)

• auto.arima() chooses AR(2) with min AICC.

• AR(2) with constant mean was fit directly to data. With observation \(Y_t\),
  \[ Y_t \hspace3mm = \hspace3mm \mu + X_t \\\\ X_t \hspace3mm = \hspace3mm \phi_1 X_{t-1} + \phi_2 X_{t-2} + e_t \]

  #- forecast ARIMA
  plot(forecast(Fit1))

  #--- 10-step prediction ---
  Y   <- Lake         # Y should be TS object
  Fit <- Fit1         # Fit you are using for forecasting
  
  Y.h  <- predict(Fit, n.ahead=10)   # Y.h$pred is TS object
  Yhat <- ts(cbind(Y.h$pred, Y.h$pred-1.96*Y.h$se, Y.h$pred+1.96*Y.h$se), start=start(Y.h$pred), freq=frequency(Y))  
  ts.plot(Y, Yhat, type="o", col=c("black","red","blue","blue"))
  abline(h=0)
  abline(h=mean(Y), col="blue")

  ##- 1. Direct ARMA Fit 

  #--- Rolling 1-step prediction ---
  Y <- Lake
  Window.size <- 60           # must be less than length(Y)
  ARIMA.order <- c(0,0,2)     # model to fit each step

  Yhat <- matrix(0, (length(Y)-Window.size), 3)         # initialize what needs to be saved 
  for (i in 1:(length(Y)-Window.size)) {
    Fit      <- Arima( Y[i:(i+Window.size-1)], order=ARIMA.order)        # <--- Fit ARMA(p,q)
    Y.h      <- predict(Fit, n.ahead=1)                                  # one step prediction
    Yhat[i,] <- c(Y.h$pred, Y.h$pred-1.96*Y.h$se, Y.h$pred+1.96*Y.h$se) 
  }
  Yhat <- ts(Yhat, start=time(Y)[Window.size+1], freq=frequency(Y))

  #- Calculate prediction performance
  Pred.error = Y[(Window.size+1):length(Y)] - Yhat[,1]
  Pred.rMSE  = sqrt(  mean( (Pred.error)^2 ) )           #- prediction root Mean Squared Error
  Pred <- matrix(c(Pred.rMSE, 1.96*sqrt(Fit$sigma2), mean(Pred.error)), 1, 3)
  colnames(Pred) <- c("Pred.rMSE", "95%PI", "Mean")
  Pred

  #- Plot the prediction result with original data ---
  ts.plot(Y, Yhat, type="o", col=c("black", "blue", "red", "red"))

  #- Plot of the prediction error with ACF of pred error
  layout(matrix(1:2, 1, 2))    # two plots side by side
  plot( Pred.error, type="o", ylim=c(-Pred[2], Pred[2])*1.2, main="Prediction Error (Blue-Black)" )
  abline(h=0)
  abline(h=c(-1.96, 1.96)*sqrt(Fit$sigma2), col="blue", lty=2)
  acf(Pred.error)
  layout(1)                    # reset the layout

##      Pred.rMSE    95%PI        Mean
## [1,] 0.8029921 1.526601 -0.07798477

2. Foreasting with ARMA with Linear Trend

If we choose that the level is going down linearly, we can try to fit a line with ARMA errors.

  Reg2 <- lm(Lake~time(Lake))
  summary(Reg2)

## 
## Call:
## lm(formula = Lake ~ time(Lake))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.50919 -0.74760 -0.01556  0.75966  2.53409 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 55.278141   7.922614   6.977 4.02e-10 ***
## time(Lake)  -0.024071   0.004119  -5.843 7.16e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.136 on 95 degrees of freedom
## Multiple R-squared:  0.2644, Adjusted R-squared:  0.2566 
## F-statistic: 34.14 on 1 and 95 DF,  p-value: 7.165e-08

  plot(Lake, type="o")                     
  abline(Reg2, col="red")         # overlay regression line in red

  plot(Reg2$residuals, type="o")            # plot of the residuals from regression

  Fit2 <- auto.arima(Reg2$residuals, d=0)   # fit the residuals with ARMA(p,q)
  Fit2

## Series: Reg2$residuals 
## ARIMA(1,0,1) with zero mean 
## 
## Coefficients:
##          ar1     ma1
##       0.6671  0.3827
## s.e.  0.0937  0.1135
## 
## sigma^2 estimated as 0.452:  log likelihood=-98.72
## AIC=203.44   AICc=203.7   BIC=211.16

  Randomness.tests(Fit2$resid)

##   B-L test H0: the sereis is uncorrelated
##   M-L test H0: the square of the sereis is uncorrelated
##   J-B test H0: the sereis came from Normal distribution
##   SD         : Standard Deviation of the series

##       BL15  BL20  BL25 ML15  ML20    JB    SD
## [1,] 0.985 0.974 0.969 0.19 0.189 0.782 0.669

  Fit2b <- auto.arima(Lake, d=0, xreg=time(Lake) )   # fit the residuals with ARMA(p,q)
  Fit2b

## Series: Lake 
## Regression with ARIMA(1,0,1) errors 
## 
## Coefficients:
##          ar1     ma1  intercept     xreg
##       0.6682  0.3817    55.5443  -0.0242
## s.e.  0.0936  0.1136    18.0324   0.0094
## 
## sigma^2 estimated as 0.4612:  log likelihood=-98.66
## AIC=207.33   AICc=207.98   BIC=220.2

Analysis 2 (linear trend)

• Linear trend was fit to \(Y_t\), then the residuals were fitted with ARMA(1,1).

• Y_t = a+b t + X_t X_t

• In other words, \[ Y_t \hspace3mm = \hspace3mm a + bt + X_t \\\\ X_t \hspace3mm = \hspace3mm \phi_1 X_{t-1} + e_t + \theta_1 e_{t-1} \]

  #- forecast ARIMA
  plot(forecast(Fit2b,  xreg=1972:1982))
  abline(a=Fit2b$coef["intercept"], b=Fit2b$coef["xreg"], col="red")

  #--- 10-step prediction ---
  Y   <- Lake         # Y should be TS object
  Fit <- Fit2b        # Fit you are using for forecasting
  
  Y.h  <- predict(Fit, n.ahead=10, newxreg=c(1:10)+max(time(Y)) )   # Y have to have freq=1
  Yhat <- ts(cbind(Y.h$pred, Y.h$pred-1.96*Y.h$se, Y.h$pred+1.96*Y.h$se), start=start(Y.h$pred), freq=frequency(Y))  
  ts.plot(Y, Yhat, type="o", col=c("black","red","blue","blue"))
  abline(h=0)
  abline(a=Fit$coef["intercept"], b=Fit$coef["xreg"], col="red")

  ##- 2. ARMA with Linear Fit 

  #--- Rolling 1-step prediction ---
  Y <- Lake
  Window.size <- 60           # must be less than length(Y)
  ARIMA.order <- c(1,0,1)     # model to fit each step

  Yhat <- matrix(0, (length(Y)-Window.size), 3)         # initialize what needs to be saved 
  for (i in 1:(length(Y)-Window.size)) {

    #- Fit ARMA(p,q) WITH Reg Line
    Fit      <- Arima( Y[i:(i+Window.size-1)], order=ARIMA.order, xreg=time(Y[i:(i+Window.size-1)]))

    #- 1-step Prediction for next Year
    Y.h      <- predict(Fit, n.ahead=1, newxreg=1+max(time(Y[i:(i+Window.size-1)])) )   
    Yhat[i,] <- c(Y.h$pred, Y.h$pred-1.96*Y.h$se, Y.h$pred+1.96*Y.h$se) 
  }
  Yhat <- ts(Yhat, start=time(Y)[Window.size+1], freq=frequency(Y))

  #- Calculate prediction performance
  Pred.error = Y[(Window.size+1):length(Y)] - Yhat[,1]
  Pred.rMSE  = sqrt(  mean( (Pred.error)^2 ) )           #- prediction root Mean Squared Error
  Pred <- matrix(c(Pred.rMSE, 1.96*sqrt(Fit$sigma2), mean(Pred.error)), 1, 3)
  colnames(Pred) <- c("Pred.rMSE", "95%PI", "Mean")
  Pred

  #- Plot the prediction result with original data ---
  ts.plot(Y, Yhat, type="o", col=c("black", "blue", "red", "red"))

  #- Plot of the prediction error with ACF of pred error
  layout(matrix(1:2, 1, 2))    # two plots side by side
  plot( Pred.error, type="o", ylim=c(-Pred[2], Pred[2])*1.2, main="Prediction Error (Blue-Black)" )
  abline(h=0)
  abline(h=c(-1.96, 1.96)*sqrt(Fit$sigma2), col="blue", lty=2)
  acf(Pred.error)
  layout(1)                    # reset the layout

##      Pred.rMSE    95%PI      Mean
## [1,] 0.7817429 1.497827 0.2225482

3. Forecasting with ARIMA (Random Trend)

Let’s look at the difference between the observations: \[ diff(Y) = Y_t - Y_{t-1} \\ \\ \]

  diff.Y <- diff(Lake)      # takes difference between the observations
  plot(diff.Y, type='o')

  layout(matrix(1:2, 1, 2))
  acf(diff.Y)
  pacf(diff.Y)

  layout(1)

  Stationarity.tests(diff.Y)       # See if diff(Y) is stationary  

## Warning in adf.test(A): p-value smaller than printed p-value

## Warning in pp.test(A): p-value smaller than printed p-value

## Warning in kpss.test(A): p-value greater than printed p-value

##        KPSS  ADF   PP
## p-val:  0.1 0.01 0.01

  Fit3 <- auto.arima(Lake, d=1, stepwise=FALSE)    # Fit the difference of Lake with ARMA 
  Fit3

## Series: Lake 
## ARIMA(1,1,2) 
## 
## Coefficients:
##          ar1      ma1      ma2
##       0.6385  -0.5349  -0.3514
## s.e.  0.1345   0.1445   0.1055
## 
## sigma^2 estimated as 0.4812:  log likelihood=-99.88
## AIC=207.76   AICc=208.2   BIC=218.02

  Randomness.tests(Fit3$residuals)

##   B-L test H0: the sereis is uncorrelated
##   M-L test H0: the square of the sereis is uncorrelated
##   J-B test H0: the sereis came from Normal distribution
##   SD         : Standard Deviation of the series

##       BL15  BL20  BL25  ML15  ML20    JB    SD
## [1,] 0.979 0.956 0.942 0.524 0.495 0.577 0.677

Analysis 3 (take difference)

• Difference of \(Y_t\) was taken. The difference seems to be WN.

• \[ \bigtriangledown Y_t \mbox{ is WN} \\ Y_t \mbox{ is ARIMA(0,1,0)} \]

  #- forecast ARIMA
  plot(forecast(Fit3))

  #--- 10-step prediction ---
  Y   <- Lake         # Y should be TS object
  Fit <- Fit3         # Fit you are using for forecasting
  
  Y.h  <- predict(Fit, n.ahead=10)   # Y.h$pred is TS object
  Yhat <- ts(cbind(Y.h$pred, Y.h$pred-1.96*Y.h$se, Y.h$pred+1.96*Y.h$se), start=start(Y.h$pred), freq=frequency(Y))  
  ts.plot(Y, Yhat, type="o", col=c("black","red","blue","blue"))
  abline(h=0)
  abline(h=mean(Y), col="blue")

  ##- 3. ARIMA Fit (Random Trend)

  #--- Rolling 1-step prediction ---
  Y <- Lake
  Window.size <- 60           # must be less than length(Y)
  ARIMA.order <- c(1,1,2)     # model to fit each step

  Yhat <- matrix(0, (length(Y)-Window.size), 3)         # initialize what needs to be saved 
  for (i in 1:(length(Y)-Window.size)) {
    Fit      <- Arima( Y[i:(i+Window.size-1)], order=ARIMA.order)        # <--- Fit ARMA(p,q)
    Y.h      <- predict(Fit, n.ahead=1)                                  # one step prediction
    Yhat[i,] <- c(Y.h$pred, Y.h$pred-1.96*Y.h$se, Y.h$pred+1.96*Y.h$se) 
  }
  Yhat <- ts(Yhat, start=time(Y)[Window.size+1], freq=frequency(Y))

  #- Calculate prediction performance
  Pred.error = Y[(Window.size+1):length(Y)] - Yhat[,1]
  Pred.rMSE  = sqrt(  mean( (Pred.error)^2 ) )           #- prediction root Mean Squared Error
  Pred <- matrix(c(Pred.rMSE, 1.96*sqrt(Fit$sigma2), mean(Pred.error)), 1, 3)
  colnames(Pred) <- c("Pred.rMSE", "95%PI", "Mean")
  Pred

  #- Plot the prediction result with original data ---
  ts.plot(Y, Yhat, type="o", col=c("black", "blue", "red", "red"))

  #- Plot of the prediction error with ACF of pred error
  layout(matrix(1:2, 1, 2))    # two plots side by side
  plot( Pred.error, type="o", ylim=c(-Pred[2], Pred[2])*1.2, main="Prediction Error (Blue-Black)" )
  abline(h=0)
  abline(h=c(-1.96, 1.96)*sqrt(Fit$sigma2), col="blue", lty=2)
  acf(Pred.error)
  layout(1)                    # reset the layout

##      Pred.rMSE    95%PI       Mean
## [1,] 0.7427918 1.498318 0.09050329

Comparing MSE

• From Model1, 2, 3, estimate for \(\sigma^2\) was \(0.4812, \, 0.4612, \, 0.4784\) respectively.

• Taking square root, that means \(\hat \sigma\) are \(0.6936, \, 0.6791, \, 0.6916\)

• On the other hand, rolling 1-step prediction rMSE for the three models were \(0.8029, \, 0.7817, \, 0.7427\)

• Why do they not agree? Which model is fitting the best?

• What is over-fittting?